
#include "cNumerov.h"

Numerov::Numerov()
{
	setName("Numerov method");
	double ivy,ivy2,ivdy, ivdy2;
	ivy=0.01;
	ivdy=0.01;
        ivy2=0.01;
	ivdy2=0.01;
	set_ty(ivy);
	set_fy(ivy);
	set_tdy(ivdy);
	set_fdy(ivdy);
        set_ivaluefy2(ivy2);
        set_ivaluety2(ivy2);
        set_ivaluefdy2(ivdy2);
        set_ivaluetdy2(ivdy2);
}

Numerov::~Numerov()
{
}

void Numerov::methodStep(double& step, vector<double> &xx, vector<double> &yy, vector<double> &dy, cDY* DY, cPot* V) {
   // std::cout << "x.back() = " << xx.back() << " y.back() = " << yy.back() << std::endl;
    yy.push_back( (( 2.0 + (5.0/6.0)  * step * step * DY->g(xx.back(), V->potential(xx.back()), yy.back()) ) * yy.back() -
				  (  1.0 - (1.0/12.0) * step * step * DY->g(xx[xx.size()-2], V->potential(xx[xx.size()-2]), yy[yy.size()-2]) ) * yy[yy.size()-2]) /
				  (  1.0 - (1.0/12.0) * step * step * DY->g(xx.back() + step, V->potential(xx.back() + step), yy.back()) )
	);
	xx.push_back(xx.back()+step);
        dy.push_back((yy.back() - yy[yy.size()-2])/step);
      //  dy.push_back((yy.back() - yy[yy.size()-2])/(2.0*step));
}

//! The method of numerical method.  This method is Numerov's method's own implementation of the "method" - method.
void Numerov::method(cDY* DY, cPot * V) {
        //Initialize all attributes
	initAttribs(V,DY);
        double nextDiscont = 0.0;
        int indx = 0;
        double step = getStep();
        
        vector<double> x,y,tx,fx,ty,fy,tdy,fdy;
        x = get_x();
        y = get_y();
        tx = getTx();
        fx = getFx();
        ty = getTy();
        fy = getFy();
        tdy = getTdy();
        fdy = getFdy();
        //cout << tx.back() << ", " << fx.back() << ", " <<  ty.back() << ", " << fy.back() << std::endl;
	//Calculate and set the gluepoint
	V->setGlue(V->calculateGluepoint(DY->get_E()));
        
        if(V->discontVectorExists()) {
            indx = V->getDiscontVector().size();
            nextDiscont = V->getDiscontVector()[indx];
        }
        else {
            nextDiscont = V->getGlue();
        }
        
        //Calculate the steps from the right side to the gluepoint.
		double nstep = -1.0*step;
                while(true) {

			try {
			    //Calculate to the gluepoint
                              if((fx.back() - V->getGlue()) < step) {
                                  if((fx.back() - V->getGlue()) > (step/100.0)) {
                                        double sstep = -1.0*(fx.back() - V->getGlue());
  //                                      std::cout << "sstep = " << sstep << std::endl;
                                        methodStep(sstep, fx, fy, fdy, DY, V);
                                  }
//                                  std::cout << "fx to gluepoint, fx.back() = " << fx.back() << std::endl;
                                break;
                            }
			    else {
			        methodStep(nstep, fx, fy, fdy, DY, V);
			    }
			} catch (cNumericalError &ne) {
				cout << "Exception occurred with message: " << endl;
				cout << ne.what() << endl;
				x.clear();
				y.clear();
				tdy.clear();
				fdy.clear();
				ty.clear();
				fy.clear();
				tx.clear();
				fx.clear();
				throw;
			}
                    indx--;
                 if(V->discontVectorExists()) {
                        nextDiscont = V->getDiscontVector()[indx];
                       }   
		}
                indx = 0;
                if(V->discontVectorExists()) {
                        nextDiscont = V->getDiscontVector()[indx];
                }
                else {
                        nextDiscont = V->getGlue();
                }
                
                while(true) {
                        try {
                              if((V->getGlue() - tx.back()) < step) {
                                  if((V->getGlue() - tx.back()) > (step/100.0)) {
                                        //Take the step to the gluepoint
                                        double tstep = (V->getGlue() - tx.back());
    //                                    std::cout << "tstep = " << tstep << std::endl;
                                        methodStep(tstep, tx, ty, tdy, DY, V);   
                                  }
                                //Take additional step over the gluepoint
                                methodStep(step, tx, ty, tdy, DY, V);   
                                //Save the additional step for later use
                                Dtx = tx.back();
                                Dty = ty.back();
                                //Remove the additional step from results
                                tx.pop_back();
                                ty.pop_back();
                                tdy.pop_back();
                                break;
                            }
			    else {
                                methodStep(step, tx, ty, tdy, DY, V);                              
			    }
			} catch (cNumericalError & ne) {
				cout << "Exception occurred with message: " << endl;
				cout << ne.what() << endl;
                                clearAllVectors();
				throw;
			}
  
                    indx++;
                     if(V->discontVectorExists()) {
                        nextDiscont = V->getDiscontVector()[indx];
                     }
		}
                
                
                setTx(tx);
                setTy(ty);
                setTdy(tdy);
                setFx(fx);
                setFy(fy);
                setFdy(fdy);
   
   //std::cout << "Before combine fy.back()-ty.back() = " << fy.back()-ty.back() << std::endl;
   combineAttribVectors();
   //std::cout << "After combine fy.back()-ty.back() = " << getFy().back()-getTy().back() << std::endl;
   
   y = normalizeInteg(get_y(),get_x());
   //std::cout << "After normalize y[ty.size()-1]-y[ty.size()] = " << y[ty.size()-1]-y[ty.size()] << std::endl;
   
   
//    double eta = fabs(ty.back()/fy.back());
    //std::cout << "ty.back() = " << ty.back() << ", " << "fy.back() = " << fy.back() << std::endl;
    //std::cout << "eta = " << eta << std::endl;
    //Dty = (eta * Dty)/getNormFactor();
    Dty = Dty/getNormFactor();
    //std::cout << "Glue is at " << tx.back() << " = " << V->getGlue() <<  std::endl;
//    
//    std::cout << "Dtx = " << Dtx << " Dty = " << Dty << std::endl;
//    
    
    //std::cout << "fabs(y[ty.size()-1] + y[ty.size()]) = " << fabs(y[ty.size()-1] + y[ty.size()]) << std::endl;
    //std::cout << "y[ty.size()-1] * y[ty.size()] = " << y[ty.size()-1] * y[ty.size()] << std::endl;
    
   if (fabs(y[ty.size()-1] + y[ty.size()]) < 0.2 && y[ty.size()-1] * y[ty.size()] < 0.0) {
             for(unsigned int indx=tx.size(); indx < y.size(); indx++){
		y[indx]= y[indx] * -1.0;
            }
	}
//   std::cout << "Dty-y[ty.size() + 1] = " << Dty << " - " << y[getTy().size() + 1] << " = " << Dty-y[getTy().size() + 1] << std::endl;
//   std::cout << "Dtx = " << Dtx << " x[getTy().size() + 1] = " << get_x()[getTy().size() + 1] << std::endl;
    setY(y);
 
     setGlueIndex(ty.size() - 1);
 
}

void Numerov::findNextEigenfunction(cDY* DY, cPot* V, unsigned break_value) {
     //double minE_Step = 1e-9;
    //TODO: Siirrettävä muualle setMinE_step();
    double D0,D1;
  //  double ee = 1e-3;
    double initialEnergy = DY->get_E();
    setMinE_step(1e-6);
    method(DY,V);
    //unsigned int bailout = 0;

    //D0 = fabs(dty/ty.back()) - fabs(dfy/fy.back());
    D0 = Dty - get_y()[getTy().size() + 1];
    DY->set_E(DY->get_E() + getE_step());
    method(DY,V);
    //D1 = fabs(dty/ty.back()) - fabs(dfy/fy.back());
    D1 = Dty - get_y()[getTy().size() + 1];

    while(fabs(D0) < fabs(D1)) {
        DY->set_E(DY->get_E() + getE_step());
        D0 = D1;
        method(DY,V);
        D1 = Dty - get_y()[getTy().size() + 1];
    }

   while(fabs(getE_step()) > getMinE_step() && DY->get_E() > initialEnergy) {

        while(fabs(D0) > fabs(D1)) {
            DY->set_E(DY->get_E() + getE_step());
            D0 = D1;
            method(DY,V);
            D1 = Dty - get_y()[getTy().size() + 1];
            //cout << "D0 = " << D0 << " D1 = " << D1 << " : fabs(D0) > fabs(D1) = " << (fabs(D0) > fabs(D1)) << endl;
        }

        setE_step(-0.5*(getE_step()));
        DY->set_E(DY->get_E() + getE_step());
        D0 = D1;
        method(DY,V);
        D1 = Dty - get_y()[getTy().size() + 1];
    }
   // cout.precision(15);
   // cout << "Energy = " << DY->get_E() << " D0 = " << D0 << " D1 = " << D1 << " tx.back() = " <<  tx.back() << " fx.back() = " << fx.back() << endl;  
}